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I. INTRODUCTION 

Systems of interacting oscillators are ubiquitous in science. 
In the common case where the natural frequencies or ampli- 
tudes or inter-oscillator couplings are time- varying, they pose 
a continuing challenge to the time-series analyst who endeav- 
ours to understand the underlying system from the signal(s) 
it creates. Oversimplifications of hypotheses are often used 
to render the problem more tractable, but can all too easily 
result in a failure to describe phenomena that are in fact of 
central importance - given that the strength, direction and 
functional relationships that define the nature of the interac- 
tions can cause qualitatively new states to appear or disap- 
pear. Time- variability of this kind is especially important in 
biological applications, though it is by no means restricted to 
biology. 

In the absence of time- variability, there are many different 
methods available 01-01 for detecting and quantifying the cou- 
plings and directionality (dominant direction of influence) be- 
tween oscillators based, especially, on the analysis of phase 
dynamics. Approaches to the detection of synchronization 
have mostly been based on the statistical properties of the 
phase difference The inference of an underlying phase 

model has been used as the functional basis for a number of 
techniques to infer the nature of the phase-resetting curves, in- 
teractions and structures of networks JH-H . However, these 
techniques inferred neither the noise dynamics nor the pa- 
rameters characterizing the noise. An additional challenge to 
these methods can be the time-varying dynamics mentioned 
above. In a separate line of development, Bayesian inference 
was applied to analyse the system dynamics lfl5r - [20l1 . thereby 
opening the door to inference of noisy time-evolving phase 
dynamics. Methods based on transfer entropy and Granger 
causality have a generality that has facilitated a number of ap- 
plications, including inference of the coupling strength and 
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directionality I21M23I1 . These techniques provide measures of 
the amount of information in a measured signal, or the causal 
relationships between measured signals and, in doing so, they 
infer effect rather than mechanism. 

In this paper we describe in detail and further extend a re- 
cently introduced method B24il based on dynamical Bayesian 
inference. As we demonstrate below, it enjoys many advan- 
tages over earlier approaches. One of these is that it does not 
require the observable to fill the domain of the probability den- 
sity function at equilibrium: it can thus provide the same in- 
formation from a small fraction of the data that is required by 
the transfer entropy or Granger causality approaches. Addi- 
tionally, the dynamical approach has the advantage of charac- 
terizing the system completely (not only in terms of informa- 
tion measures). Thus, from the inferred dynamics one can 
deduce self-consistently any information that is of interest, 
be it coupling functions, or synchronization, or causality, or 
equilibrium densities, including the equations of motion. We 
discuss in detail the theoretical background, the technical as- 
pects and limitations of the algorithms, and we demonstrate 
the wide applicability of the method by consideration of sev- 
eral examples. 

The coupling functions are of particular importance. Their 
form is uniquely able to describe the functional laws of in- 
teraction between the oscillators. Earlier theoretical studies 
have included the work of Kuramoto |25] and Winfree i26lL 
which used a function defined either by the phase difference 
or by both phases, and of Daido Hzl Q and Crawford fH 
who used a more general form in which the coupling func- 
tion was expanded in Fourier series. Other methods for in- 
ference of the coupling functions have also been suggested 
EStmmiliQl]- The technique described below goes beyond 
all of these because it is able to follow the time- variability of 
the coupling functions and hence can reveal their dynamical 
character where it exists. 

We will also show how the technique can readily be 
extended to encompass networks of interacting oscillators. 
These form a large and important grou p of physical systems, 
including neural networks J11 l3ll l32ll. electrochemical sys- 
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terns ifToL l33ll , crowd synchrony on the Millennium bridge, 
and networks of fireflies [3 . The large scale of the networks 
can introduce a higher complexity, both in structure and func- 
tional behavior. For example in neuronal networks, the ex- 
istence of spatial and spatial-temporal correlations, collective 
or partially collective (clustering) behavior, synchronization 
or desynchronization, and time- variability has been reported 
Oil |32[ |35[ |36|]. In such cases, and given the kinds of phe- 
nomenon to be studied, there is an increasing need for power- 
ful techniques that can infer the time- varying dynamics of the 
oscillatory networks. 

In Sec. HJ we provide details about the phase decomposi- 
tions, the implementation of the Bayesian framework and how 
the time- varying information is propagated. The synchroniza- 
tion detection through a map representation of the phase dy- 
namics is discussed in Sec. Hill while the method for describ- 
ing the interactions is demonstrated in Sec. |lVl Before the 
method is applied, we consider in Sec.[V]some important tech- 
nical aspects and limitations. The wide applicability of the 
method is demonstrated in Sec. [VH through the analysis of 
time-series from numerical phase and limit-cycle oscillators, 
analogue simulation and cardio-respiratory interactions. The 
generalization of the approach to networks of oscillators, as 
exemplified by two numerical examples, is presented in Sec. 
VIII Finally, we summarise and draw conclusions in Sec. IVIII1 
The algorithm used for the detection of synchronization is de- 
scribed in Appendix lAl 



n. PHASE-DYNAMICS DECOMPOSITION 

Consider an TV-dimensional oscillator dx/dt = f(x(t)) 
whose solution f admits a limit cycle. Such an oscillator can 
usually be represented by a constant phase velocity <j> = 00 and 
a vector coordinate that defines the limit cycle as a function of 
the phase </>: r = r ((/>). 

When two such oscillators mutually interact sufficiently 
weakly, their motion is commonly approximated just by their 
phase dynamics 111, [37|]. We note that, in general, if we 
describe the phase of a system through a generic monotonic 
change of variables, than the dynamical process can be writ- 
ten as 



<t>i=UJi + fi((/>i) +j 



(1) 



Eq. (Q]) explicitly includes a noise term & to enable it to repre- 
sent a process in a real system. The noise can be e.g. a white 
Gaussian noise (£i(t)£j(r)) = S(t — r)Eij , which may pos- 
sibly be spatially correlated. 

The phase-dynamics decomposition technique is highly 
modular from the algorithmic point of view, and each mod- 
ule will be explained separately in the sections that follow. 
The overall procedure comprises the following steps: 

1. Assumption that the dynamics can be precisely de- 
scribed by a finite number of Fourier terms (see Sec. 
UTAb. 



probability distribution (see Sec. Ill B I for stationary dy- 



namics, and Sec. Ill CI for time- varying dynamics). 

3. Integration of the probability that this parameter set lies 
inside the Arnold tongue defining synchronization. This 
effectively yields the cumulative probability of the syn- 
chronization state of the dynamics (see Sec. IIII Al) . 

4. Use of the parameter information as obtained in step 2 
to create a description of the interactions, leading to de- 
tection of the predominant directionality and coupling 
function estimation among the oscillators (see Sec. HVl) . 



A. Truncated Fourier Series 

The periodic behaviour of the system suggests that it can 
appropriately be described by a Fourier decomposition. De- 
composing both fi and gi in this way leads to the infinite sums 
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and 



00 00 



i(<f>i, 4>j) = Yl E c i;r , s e a ^e i2 ^. (2) 
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It is reasonable to assume that, in most cases, the dynamics 
will be well-described by a finite number K of Fourier terms, 
so that we can rewrite the phase dynamics of Eq.(Q} as a finite 
sum of base functions 

K 

4>\ = + (3) 

k=-K 

where 1=1,2, $1,0 = $2,0 = 1, $ — and the rest of 
and are the K most important Fourier components. 



B. Bayesian Inference 

In order to reconstruct the parameters of Eq. ([3]) we make 
extensive use of the approach to dynamical inference pre- 
sented in (H, [l^] . In this section we briefly outline the tech- 
nique as adapted to the present case. The fundamental prob- 
lem in dynamical inference can be defined as follows. A 2- 
dimensional (in general L-dimensional) time-series of obser- 
vational data X = {4>i, n = &(tn)} (t n = nh) is provided, 
and the unknown model parameters M = {c%\ , } are to 
be inferred. 

Bayesian statistics employs a given prior density p pr ior(A1) 
that encloses expert knowledge of the unknown parameters, 
together with a likelihood function t(X\M), the probability 
density to observe {4>i, n (t)} given the choice M of the dy- 
namical model. B ayes' theorem 



2. Inference, given the data, of the Fourier terms, the noise 
amplitude, and their correlation in form of a parameter 



Px(M\X) 



£(X\M) Ppnor (M) 

f£(X\M) Pv nor(M)dM 



(4) 
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then enables calculation of the so-called posterior density 
Px(-M\X) of the unknown parameters M conditioned on the 
observations. 

For independent white Gaussian noise sources, and in the 
mid-point approximation where 



»l,n 



0Z,n+l — <l>l,n 

h 



and 



k,n + 0Z,n+l) 



,J l,n 



the likelihood is given by a product over n of the probabil- 
ity of observing 0z, n +i at each time. The likelihood function 
is constructed by evaluation of the stochastic integral of the 
noise term over time, as 



&\U)= [ t%+1 £ l (t)dt=VhH 
Ju 



(5) 



where H is the Cholesky decomposition of the positive defi- 
nite matrix E, and z\ is a spatially uncorrelated standardized 
Gaussian variable. The joint probability density of z\ is used 
to find the joint probability density of the process in respect of 

(MU+i)-MU))byim V osmgP(fa(t i+1 ) = det(J^)P(f), 

where is the Jacobian term of the transformation of vari- 
ables that can be calculated from Eq. (|2]). If the sampling fre- 
quency is high enough, the time step h tends to zero, and the 
determinant of the Jacobian can be well- approximated by 
the product of its diagonal terms 



n 



d<$> 



iM 



This transformation leads to an extra term in a least squares 
likelihood, and the minus log-likelihood function S = 
— hii(X\M) can thus be written as 



N 



S=-\n\E\ + - £(. fc 

n=0 



(6) 



where summation over the repeated indices k,l,i ,j is implicit. 

The log-likelihood © is a quadratic form of the Fourier 
coefficients of the phases. Hence if a multivariate prior prob- 
ability is assumed, the posterior probability is a multivariate 
normal distribution as well. 

This is highly desirable for two reasons: (i) a Gaussian pos- 
terior is computationally convenient because it guarantees a 
unique maximum, with the mean vector and covariance ma- 
trix completely characterizing the distribution and giving us 
the most significant information; (ii) all the multivariate nor- 
mal posteriors can be used again as priors in the presence of a 
new block of data, and knowledge about the system can eas- 
ily be updated. This last feature is essential for any real-time 
application because it ensures that the complexity of the algo- 
rithm does not change with the length of the input data-stream. 

From and assuming a multivariate normal distribution 
as the prior for parameters , with means c, and co variances 



Ep r i or , the stationary point of S can be calculated recursively 
from 
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c (i) = fc-l-vW) r (l) 

2 d<j>i ' 



(7) 



where the covariance is XI = H _1 , summation over n from 
1 to TV is assumed and summation over repeated indices 
k,l,i,j,wis again implicit. 

We note that a noninformative "flat" prior can be used as 
the limit of an infinitely large normal distribution, by setting 

"prior = and Cpnor = 0. 

The multivariate probability Mx (c|, c, 5) given the readout 
time series X = {0z, n = ^i(t n )} explicitly defines the prob- 
ability density of each parameter set of the dynamical system. 
Because each of them can be discriminated, as belonging or 
not belonging to the Arnold tongue region we can define the 
binary property s(cj^) = {1,0}, and can obtain the posterior 
probability of the system being synchronized or not by evalu- 
ating the probability of s 

Psync=Px(s=l) = J s(c) M X (c|c, S) dc . (8) 

The computation of p sync will be discussed in Sec. [TTTJ 

C. Time- varying information propagation 

The multivariate probability described by Nx{c, S) for the 
given time series X = {(j) n = (j>(t n )} explicitly defines the 
probability density of each parameter set of the dynamical 
system. When the sequential data come from a stream of 
measurements providing multiple blocks of information, one 
applies © to each block. Within the Bayesian theorem, the 
evaluation of the current distribution relies on the evaluation 
of the previous block of data, i.e. the current prior depends on 
the previous posterior. Thus the inference defined in this way 
is not a simple windowing, but each stationary posterior de- 
pends on the history of the evaluations from previous blocks 
of data. 

In classical Bayesian inference, if the system is known to 
be non-time- varying, then the posterior density of each block 
is taken as the prior of the next one: = Sp 0St . This full 

propagation of the covariance matrix allows good separation 
of the noise, and the uncertainties in the parameters steadily 
decrease with time as more data are included. 

If time- variability exists, however, this propagation will act 
as a strong constraint on the inference, which will then fail 
to follow the variations of the parameters. This situation is 
illustrated in Fig. [Ha) ll53ll . In such cases, one can consider 



the processes between each block of data to be independent 
(i.e. Markovian). There cannot then be any information prop- 
agation between the blocks of data, and each inference starts 
from the flat distribution = 0. The inference can thus 

follow more closely the time- variability of the parameters, but 
the effect of noise and the uncertainty of the inference will of 
course be much larger, as shown in Fig. [Tib). 

Where the system's parameters are time-dependent, we 
may assume that their probability diffuses normally accord- 
ingly to the known diffusion matrix E^iff- Thus, the proba- 
bility density of the parameters is the convolution of the two 
normal multivariate distributions, E post and Xldiff 

v n +i V n _i_ V n 

^prior — ^post "•" ^diff' 

The covariance matrix Ediff expresses our belief about which 
part of the dynamical fields that define the oscillators has 
changed, and the extent of that change. Its elements are 
(^diff)i,j = Pij^iCJj, where is the standard deviation of the 
diffusion of the parameter q after the time window t w that has 
elapsed from the first block of information to the following 
one. pij is the correlation between the change of the parame- 
ters Ci and Cj (with pa = 1). In relation to the latter, a special 
example of Ediff will be considered: we assume that there is 
no correlation between parameters, i.e. pij =0, and that each 
standard deviation di is a known fraction of the parameter q : 
&i = PwCi (where p w indicates that p is referred to a window 
of length t w ). It is important to note that this particular exam- 
ple is actually rather general because it assumes that all of the 
parameters (from the Ep 0St diagonal) can be of a time- varying 
nature - which corresponds to the inference of real (experi- 
mental) systems with a priori unknown time- variability. 

There are two obvious limits in modeling the knowledge 
assumed with respect to possible time variation of param- 
eters. The first of these is to assume no time- variability: 
in this case the full information propagation matrix is used, 

Spdor = Impost- If me assumption proves wrong, the in- 
ferred parameters may accumulate a bias when the real sys- 
tem varies in time. The other limit is to assume each time 
window to be completely independent of the previous signal 
history. In this case no propagation is used, Xlpdor = 00 > 
(i.e. Sprior = 0)' anc * there is no bias but, because much infor- 
mation is forgotten, the probability of the inferred parameters 
has a large covariance matrix. An optimal assumption must 

lies in between these two limits: £prtor = Epost + Ediff- 
If a diffusion matrix is assumed, we allow the method some 
freedom for the time- variability to be followed, while restrict- 
ing it to be unbiased. The amount of variability is part of 
the model, like the number of free parameters in any standard 
method. Fig. [T] illustrates the two extreme limits, and a possi- 
ble trade-off. The inference in Fig. [He) demonstrates that the 
time- variability is captured correctly and that the uncertainty 
is reduced because more data have been included. 

If one knows beforehand that only one parameter is vary- 
ing (or, at most, a small number of parameters), then Ediff can 
be customized to allow tracking of the time- variability specif- 
ically of that parameter. This selective propagation can be 
achieved if, for example, not all but only the selected correla- 
tion pa from the diagonal has a non-zero value. 
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Figure 1: Inference of a rapidly time-varying coupling parameter 
from coupled noisy oscillators i\2\ . The gray lines represent the ac- 
tual parameter in the numerical simulation, whereas the black lines 
indicate the time-varying parameter inferred from the resultant time 
series, for: (a) full propagation, Sp^ 1 = Sp 0St ; (b) no propaga- 
tion, Sp^ 1 = oo; and (c) propagation for time-varying processes, 

^prior — ^post i ^diff • 



III. SYNCHRONIZATION DETECTION 

It is important to note that finite noise can induce phase slips 
in a system that would be synchronized in the noiseless limit. 
Rather than focusing on the presence and statistics of phase- 
slips, we propose to detect synchronization from the nature 
of the phase- slip itself. A novel feature of the present study is 
that it proposes evaluation of the probability that the equations 
driving the dynamics are intrinsically synchronized and thus 
of whether any phase-slips that may possibly be observed are 
dynamics-related or noise-induced. 

After performing the inference, one can use the recon- 
structed parameters, derived in the form of a multivariate 
normal distribution Afx(c, £), to study the interactions be- 
tween the oscillators under study. In general, the border of 
the Arnold tongue may not have an analytic form and, even if 
it has, the integral may not have an analytic solution but must 
be evaluated numerically. In practice, we estimate p sync nu- 
merically, sampling from the parameter space many realiza- 
tions {c^} m , where m labels each parameter vector tested. 
For every set of c we compute s(c m ) numerically. Let us as- 
sume for now that s(c m ) is given. To find p sync with arbi- 
trary precision, it is enough to generate a number M of pa- 
rameters c m = {c^}m, with m = 1, . . . , M sampled from 

Nx{c\c, E), since p sy nc lim M ^oo E! s ( c m)- 

However, this 2 if -dimensional integration quickly be- 
comes inefficient with an increasing number of Fourier com- 
ponents. Moreover, as we will discuss in Sec. lIII Al the com- 
putation time of the variable s(c m ) is not insignificant. On the 
other hand, if the posterior probability px is sharply peaked 
around the mean value c, then p sync will be indistinguishable 
from s(c), and the evaluation of s(c) will suffice. 



A. Synchronization Discrimination and map representation 

We now illustrate a simple numerical technique to recog- 
nize whether a coupled phase oscillator system is synchro- 
nized, or not. The technique itself amounts to a simple check 
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Figure 2: Torus representation of the phase dynamics, with toroidal 
coordinate C(0i(*)> fait)) an ^ polar coordinate ip(<fii(t), fait)). 
The white circle denotes the Poincare cross section. 



by numerical integration of the system of ordinary differen- 
tial equation defined by Eq. (Q} through one cycle of the dy- 
namics, and testing whether the 1 : 1 synchronization condition 
\ip(t)\ = \ fait) - fait)\ < K is always obeyed. 

Let us assume we are observing motion on the torus T 2 
defined by the toroidal coordinate (ifa (£), fa it)) = ifa it) + 
fa it))/ 2, and the polar coordinate ijj(t). 

For assessment of possible 1 : 1 synchronization the phase 
difference ip(t) will be defined as ipifait), fait)) = fait) — 
fait). Fig. [2 provides a schematic representation of the phase 
dynamics on the torus. Let us consider a Poincare section de- 
fined by ( = and assume that ((t)/dt\^=o > for any ip. 
This means that the direction of motion along the toroidal co- 
ordinate is the same for every point of the section. Ideally we 
would follow the time-evolution of every point and establish 
whether or not there is a periodic orbit; if there is one, and if 
its winding number is zero, then the system is synchronized. 
If such a periodic orbit exists, then there is at least one other 
periodic orbit, with one of them being stable and the other 
unstable. 

The solution of the dynamical system over the torus yields 
a map M : [0, 2tt] —> [0, 2tt] that defines, for each ijj n on the 
Poincare section, the next phase ^n+i after one circuit of the 
toroidal coordinate ^n+i = Miip n ). Fig. Ob), (c) illustrates 
the map M as evaluated computationally in two situations, 
corresponding to no synchronization, or synchronization, re- 
spectively. 

The map M is continuous, periodic, and has two fixed 
points (one stable and one unstable) if and only if there is a 
pair of periodic orbits for the dynamical system, i.e. synchro- 
nization is verified if i(j e exists such that i(j e = M(^ e ) and 

' dM(ip) I 

dip IV'e 



< 1. 



The existence of the fixed point ijj e is estab- 
lished through the simple algorithmic procedure described in 
Appendix lAl 



IV. DESCRIPTION OF THE INTERACTIONS 

Inferring the parameters of the system not only allows for 
evaluation of the synchronization as an epiphenomenon in its 
own right, but their probability A/#(c, E) also describes the 
interaction properties of the oscillators. Because their dynam- 
ics is reconstructed separately, as described by Eq. (Q]), use 
can be made only of those inferred parameters from the base 



functions fiifij) and giifa^fa) that are linked to influences 
between the oscillators. 

One can seek to determine the properties that characterize 
the interaction in terms of a strength of coupling, predomi- 
nant direction of coupling, or even by inference of a coupling 
function. The analysis of information propagation allows in- 
ference of the time-varying dynamics, and the interactions' 
properties can be traced in time as well. This is especially im- 
portant for the inference of open interacting oscillatory pro- 
cesses where the time- variability of the interactions can lead 
to transitions between qualitatively different states, such as 
synchronization ll37ll or oscillating death ll38ll . 

The coupling amplitude quantifies the total influence be- 
tween the oscillators in a particular direction, e.g. how much 
the dynamics of the first oscillator affects the dynamical be- 
havior of the second oscillator (1 —> 2). Depending on 
whether the coupling is in only one direction, or in both di- 
rections, we speak of unidirectional or bidirectional coupling, 
respectively. In the inferential framework that we propose, the 
coupling amplitudes are evaluated as normalized measures, 
based on the interacting parameters inferred from the coupling 
base functions qiifa^fa). The quantification is calculated as 
a Euclidian norm: 



£21 
ei2 



mm, fa) 

\Q2ifa,fa) 



(9) 



where the odd inferred parameters are assigned to the base 
functions qiifa,fa) for the coupling that the second oscillator 
imposes on the first (e2i : 2 — >> 1), and vice versa (ei2 : 1 —> 
2). 

The directionality of coupling ^ often provides useful in- 
formation about the interactions. It is defined as normalization 
about the predominant coupling amplitude 



D 



€12 — e 21 
Cl2 + ^21 



(10) 



If D G (0, 1] the first oscillator drives the second (1 —> 2), or 
if D G [—1,0) the second (2 — » 1) drives the first. The quan- 
tified values of the coupling strengths or the directionality 
D represent measures of the combined relationships between 
the oscillators. Thus, a non-zero value can be inferred even 
when there is no interaction. Such discrepancies can be over- 
come by careful surrogate testing (HI S3] - by rejection of 
values below an surrogate acceptance threshold, which can be 
specified e.g. as the mean plus two standard deviations among 
many realization of the measure. 

In addition to the coupling strength and the directionality, 
one can also infer the coupling function that characterizes the 
interactions, i.e. the law that describes the functional relation- 
ships between the oscillators. Its characteristic form reflects 
the nature of the oscillators and how their dynamics reacts to 
perturbations. 

The coupling function should be 27r-periodic. In the in- 
ferential framework under study, the coupling functions were 
decomposed into a finite number of Fourier components. 
The function describing the interactions between the two 
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Figure 3: (Color online) Synchronization discrimination for the coupled phase oscillators dTTT) . (a) Schematic of an Arnold tongue in the 
coupling-frequency e-uo plane: synchronization exists only within the shaded area H3VTI . (b) Map of M(ip) for e±2 — 0.25 demonstrating that 
the oscillators are not synchronized, (c) Map of M(ip) for a case where a root of M(ip) = ip exists, i.e. where that the state is synchronized, 
(d) The corresponding phase difference, exhibiting two phase slips. 



oscillators was decomposed by use of the odd parameters 
#i(0i? 02) £ {ci, C3, . . .} and the corresponding base func- 
tions $ n [gi(0i, 2 )] e {sin(0i, 2 ), cos(0i, 2 )} up to order 
n of the decomposition. The reverse function g 2 (0i,0 2 ) £ 
{c 2 , C4, . . .} was similarly decomposed. 



V. TECHNICAL ASPECTS AND CONSIDERATIONS 

The technique is quite generally applicable to a broad class 
of problems D4IL and so there are a number of technical as- 
pects and choices to bear in mind. We now discuss three of 
them in particular: the number of base functions to be em- 
ployed in the inference process (which is part of the model); 
the intensity of the noise characterizing the system (which is 
an externally imposed constraint); and the time resolution. 

a. Number of base functions. Selection of the optimal 
set of base functions to describe the problem is far from triv- 
ial. In general one wishes to have the minimal set that de- 
scribes sufficiently well the model to be tested. Where the 
length of the data series is very long or effectively infinite, 
one can include an excessive number of base functions with- 
out immediate penalties. In reality, however, any unneeded 
base function jeopardizes the precision of the coefficients that 
really are relevant for the model, and the picture is further 
complicated when the model to be adopted is expected to be 
an outcome of the inference machinery. Where one deals with 
a long data series, possibly with a high signal-to-noise ratio, 
a relative large number of base functions can be used. T^ 
speed of computation is also an important aspect to keep un- 
der consideration, given that having a large number of base 
functions vastly increases the parameter space, and that itera- 
tive calculations (especially matrix inversion) slow the speed 
of processing by the third power of the number of coefficients. 
Note that, even though the Bayesian inference is generally 
popular in real-time applications, computational speed limi- 
tations mean that our inference framework for general phase 
dynamics cannot yet be used in this way. 

b. Role of noise intensity. In general, the greater the 
noise intensity, the bigger the covariance of the inferred pa- 
rameters. For a repeated experiment (e.g. generation of a syn- 



thetic signal, and parameter inference based on that signal) 
the variance of a particular parameter would increase mono- 
tonically with noise amplitude, as shown in Fig. [4] There are, 
however, a few notable exceptions. The inferential capabili- 
ties rely on the volume of phase- space spanned by the vari- 
ables. A state of synchronization would represent a limit cy- 
cle for the global system, and parameter inference of neither 
oscillator would reach satisfactory precision. In such cases 
a minimal amount of noise is typically needed, sufficient to 
drive the system out of equilibrium at least once. During the 
resultant phase- slip, the data would be filling the phase space 
sufficiently for correct parameter reconstruction. 

c. Time resolution. We now summarize the limits of an 
idealized data acquisition. The time step h is much bigger 
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Figure 4: Statistics of the inferred frequency uo\ and coupling e\ for 
different noise intensities E. The signal to be analysed is generated 
from Eq. (fT2l . The dotted line shows the actual values of the pa- 
rameters. The boxplots refer to the descriptive statistics (median, 
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Figure 5: Inference of (a) a time- varying frequency and (b) coupling 
parameter from time series data generated by model fi~2t for four 
different lengths of the inference windows. The sizes of the windows 
are shown in the box. 



than the correlation time of the noise (assuming the latter to 
be non-zero). At the same time, however, h is much smaller 
than any of the sequential time- windows used as data blocks 
for inference, so that each block contains many data points. 
Also, h is much smaller than either of the oscillator periods. 
Each inference block is big enough to contain many cycles 
of the dynamics (in particular, more cycles that those typical 
of a phase slip) while, at the same time, each block is small 
enough to provide the desired resolution of parameter change. 

It can happen that the time resolution of the change in dy- 
namical parameters is incompatible with an acquisition time- 
window that would guarantee precision for other parameters. 
The choice of the time- window must therefore be done on a 
case-by-case basis, depending on the type of information that 
is required from the system. Fig. \5\ illustrates such a com- 
promise. We use the numerical model (TT2t that will be intro- 
duced in Sec. I VI A"2l to investigate the time-resolution for the 
case where the frequency uj\ (t) = uj\ + A\ sin(cDt) and cou- 
pling amplitude Ei($) = £2 + A\ sin(u)t) were varying peri- 
odically at the same time. The parameters were: uj\ = 2n 1.1, 
uj2 = 2tt2.77, £1 = 0, £ 2 = 1, Cj = 2tt 0.002, A x = 0.1 
A 2 =0.5 and noise strengths E\ = E 2 = 0.15 . The parame- 
ters were reconstructed using four different window lengths 
for the inference. The results presented in Fig. \5\ demon- 
strate that, for small windows (0.5s), the parameters are sparse 
and sporadic, while for very large windows (100s) the time- 
variability is faster than the size of the window and there is 
cut-off on the form of the variability. The optimal window 
length will lie between these two. Another interesting feature 
is that, for the smallest window (0.5 s), the coupling amplitude 
improves with information propagation as time progresses, 
while the frequency inferred (as a constant component with- 
out base function) remains sparse throughout the whole time 
interval. 



VI. APPLICATIONS 

The technique is first applied to synthetic data to test the 
performance of the algorithm, and then real data are analyzed. 
To create the synthetic data, we used both numerical and ana- 
logue electronic simulations. 



A. Numerically-generated test data 

Numerically generated data were obtained from models of 
phase oscillators and limit-cycle oscillators. 



1. Phase oscillators 

The phase oscillator model provides a sufficient basis for 
the description of synchronization while being, at the same 
time, analytically traceable. We thus test the detection of 
synchronization (as explained in Sec. [TTT1) through Bayesian 
inference of synthetic data whose synchronization is already 
known. The model is given by two coupled phase oscillators 
subject to white noise 

<f>i=(jJi-\- 6jism(^j - <j>i) + &(t) ij = 1,2. (11) 

The parameters are uj\ = 1.2, 002 = 0.8, 621 =0.1; param- 
eter 612 is chosen so that the system lies close to the border 
of the Arnold tongue (either just inside or just outside). Be- 
cause we aim to demonstrate the precision of synchroniza- 
tion detection, we add no time-variability to the model, the 
inference is applied to a single block of data, and there is 
no spatial noise correlation with noise intensities En = 2. 
The dynamics of the phase difference is described as ip = 
Acj — esin^) + + &(t), where Auj = U2 — w\ is the 
frequency mismatch and e = 621 + €12 is the resultant cou- 
pling. It is evident that the analytic condition for synchroniza- 
tion, i.e. the existence of a stable equilibrium solution ip < 0, 
is Acj/e < 1. For ei2 = 0.25 (outside the Arnold tongue) 
the reconstructed map M(ip) (Fig. [3b)) after parameter in- 
ference has no root M(ip e ) = ip e : hence the oscillators are 
not synchronized. When €12 = 0.35, even though the sys- 
tem was inside the Arnold tongue, noise triggered occasional 
phase slips (see Fig.Od)). We tested synchronization detec- 
tion on the same signals using the methods already available 
in the literature, based on the statistics of the phase difference 
but none of them was able to detect the presence of 
synchronization under these conditions. In spite of the phase 
slips, our technique correctly detects the root M(ip e ) = ip e 
from the inferred parameters, revealing that the oscillators are 
intrinsically synchronized as shown in Fig. Efc): the phase 
slips are attributable purely to noise (whose inferred intensity 
is given by the matrix Eij), and not to deterministic interac- 
tions between the oscillators. 
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Figure 6: (Color online) Extraction of time- varying parameters, syn- 
chronization and coupling functions from numerical data created by 
<TT2l) . The frequency cji(t) (a) and coupling S2(t) (b) are indepen- 
dently varied. The dotted and full lines plot the parameters when the 
two oscillators are synchronized for part of the time {s\ = 0.3), and 
not synchronized at all (si = 0.1), respectively. The regions of syn- 
chronization, found by calculation of the synchronization index, are 
indicated by the gray shaded regions, (c) and (d) show the coupling 
functions q± (0i, fa) and #2(^1, ^2) for time windows centered at 
t = 350s. In both cases, the window length was t w = 50s and the 
coupling was £12 = 0.1. 



2. Limit-cycle oscillators 

To demonstrate the capabilities of the technique in trac- 
ing time- varying parameters, coupling functions, directional- 
ity and synchronization, we analyzed data from a numerical 
model of two coupled, non-autonomous, Poincare oscillators 
subject to white noise, 

±i -nxi -Ui(t) yi +6i(t) qi(xi,Xj,i) +&(£), 

y% = -nyi + Xi +e i (t)q i {y i ,y j ,t) + q 2 ) 

ri = (y/a% + y?-l) i,j = l,2. 

We tested several possibilities for the parameters: while let- 
ting the frequencies uj{ and coupling parameters Si be time- 
varying, we ran numerical experiments with the coupling 
function qi either fixed or time- varying. 

As a first numerical experiment, we considered bidirec- 
tional coupling (1^2), where the natural frequency of the 
first oscillator, and its coupling strength to the second one, 
vary periodically at the same time: ui(t) = uj\ + A\ sin(u)it) 
and £2(t) = £2 + A2 sm.{Co 2 t). The other parameters were: 
£2 0.1, uji = 2tt1, uj 2 27T1.14, A 1 = 0.2, A 2 0.13, 
uji = 2tt 0.002, Co 2 = 2tt 0.0014 and noise E u = E 22 = 0.1. 
The coupling function was proportional to the difference in 
the state variables: qi(xi,Xj,i) = — Xj and qiiyi, Vj,t) = 
Vi — Vj (the same coupling function was used for construction 
of Fig. [T] Fig. [4] and Fig. [5]). The phases were estimated as the 
angle variable fa = arctan(^/xi). With £\ = 0.1, in a state 
of no synchronization, the time- varying parameters u)\(t) and 
£ 2 (t) are accurately traced as it can be seen in Fig.fSJa) and 
(b). The coupling amplitude of £\ = 0.3 corresponds to a 
state of intermittent synchronization, where the two oscilla- 
tors are synchronized for part of the time. The precision of 



the reconstructed time- variable parameters is satisfactory dur- 
ing the non- synchronized intervals. During the synchronized 
intervals, however, the oscillators do not span sufficient phase- 
space to allow precise inference of the parameters (Fig. Oa) 
and (b), dashed lines). Within these synchronized intervals, 
the posterior probability distribution of the parameters was not 
peaked; however, it was sensibly different from zero only in 
that parameter region for which the corresponding noiseless 
dynamics is synchronized. Hence, despite the impossibility 
of accurate parameter tracking, the detection of a synchro- 
nized state (s(c) = 1) is always precise (Fig. Oa) and (b), 
grey shaded regions). 

The reconstructed sine-like functions qi{(j)i,(j) 2 ) and 
#2(01? ^2) are shown in Figs. Oc) and (d) for the first and 
second oscillators, respectively. They describe the functional 
form of the interactions between the two Poincare systems in 
Eq. JT2l) . The reconstructed form of the coupling functions 
was evaluated dynamically for each block. 

Next, the method was applied to deduce the predominant 
direction of coupling as specified from the norm of the in- 
ferred coupling base parameters. To illustrate the precision 
of the directionality detection, the frequencies were now set 
constant, while both of the coupling strengths remained dis- 
cretely time- varying. The parameters were oj\ = 27rl.3, 
UJ2 = 27T 1.7, En = ^22 =0.2, and the coupling functions 
were, as in the previous example, qi(xi,Xj,t) = X{ — Xj and 
Qi (Ui iVh^) ~ Vi~Vr Synchronization was not reached, how- 
ever, for these parameters. The couplings undergo changes at 
particular times, but otherwise remain constant, as shown in 
Fig. [71 The detected directionality index D was consistent 
with the actual values. Note that, for unidirectional coupling, 
D does not quite reach unity on account of the noise. 

To further investigate the ability to track subtle changes of 
time- varying coupling functions, we used the same model as 
in Eq. (IT2l) to generate a synthetic signal where the coupling 
functions are absolute values of the state difference to a power 
of the time- varying parameter: 

Qx,i( x i, Xj,t) = \{Xj — Xi) ()|, 

q y ,i(yi,yj,t) = \(yj -yiY {t) \, 
where z,j = {1,2} and i ^ j. The exponent parameter 
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Figure 7: Directionality of coupling D for discretely time- varying 
coupling amplitudes s\ and £2. Different unidirectionally and bidi- 
rectionally coupled states are reached for different values of s\ and 
£2, as indicated by the square insets. 
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Figure 8: (Color online) Time-evolution of the coupling function from model (fT2l) with exponential time variations JT3l ). (a)-(d) Coupling 
function qi(<pi, $2) from the first oscillator for four consecutive time windows. The window length was t w = 50s. For simplicity and clarity 
only the function q\ (<f>\ , (^2) is shown. The behavior of q2 (<pi , $2) from the second oscillator was similar. 



varied linearly with time v(t) = {1 —> 3}, and the other 
parameters remained constant: uj\ = 2n 1, uj 2 = 27r2.14, 
ei = 0.2, e 2 = 0.3 and E n = E 22 = 0.05. The recon- 
structed phase coupling functions 02 ) were calculated 
from the inferred parameters for the interacting terms of the 
base functions. The results for four consecutive windows are 
presented in Fig. [8] It can readily be seen that their complex 
form now is not constant, but varies with time. Comparing 
them in neighboring (consecutive) pairs: (a) and (b), then (b) 
and (c), then (c) and (d), one can follow the time-evolution of 
the functional form. Even though we can follow their time- 
variability, the two most distant functions Fig. [St a) and (d) are 
of substantially different shapes. Note also that, beside their 
form, the functions' norm i.e. coupling strength also varies 
(cf. the height of the maxima in Fig.^a) and (d)). 



were e 0.7, ui 2tt 15.9, ou 2 2tt 17.5, A 1 0.03, oj 
2tt 0.2 and c = 100 is constant resulting from the analogue in- 
tegration. The phases were estimated as fa = arctan(±i/;ri). 

For the given parameters the oscillators were synchronized, 
so that the second driven oscillator changed its frequency from 
being constant to being time- varying. Applying the inferential 
technique showed, correctly, that the oscillators were indeed 
synchronized (s(c) = 1) throughout the whole time period. 
The frequency of the driven oscillator was inferred as being 
time- varying Fig.[9lc). Performing a simple FFT (Fig. [3d)) 
showed that uj 2 (t) is periodic with period 0.2 Hz (exactly as 
set on the signal generator). Clearly, the technique reveals 
information about the nature and the dynamics of the time- 
variability of the parameters - and is able to do so using a 
more realistic signal than that from a numerical simulation. 



B. Analogue simulations 

We also tested the technique on signals emanating from 
analog models. These are real, highly controllable, oscillatory 
systems and the noise on their signals is real rather than con- 
trived, as in the case of numerical models. It is attributable 
to environmental disturbances, thermal fluctuations, and the 
inherent nonidealities of the circuit components. During the 
process of data acquisition and discretization, measurement 
noise can be introduced as well - noise which has no links 
with the actual dynamics of the interacting oscillators. Such 
signals provide a good test of our analysis capabilities. 

We analyzed data from an analog experimental simulation 
of two coupled van der Pol oscillators. Details of the elec- 
tronic implementation are given elsewhere j4lh . The noise 
here arises mainly from the imperfections of the electronic 
components and there is also measurement noise. 

Fig.Oa) shows the phase portrait derived from the first os- 
cillator, with time- varying frequency, which drives the second 
oscillator 

±xi - /ii(l - x\)\x x + [cji +£j 1 (t)} 2 x 1 = 0, 
4^2 - ^2(1 - xl)\x 2 + u\x 2 + e{x x - x 2 ) = 0, (14) 

where the periodic time- variability &\(t) = A\ sin(a)t) (Fig. 
[9jt>)) comes from an external signal generator. The parameters 



C. Cardiorespiratory interactions 

Having tested our technique on two quite different kinds of 
synthetic data, we now apply it to a real physiological prob- 
lem, to investigate the cardiorespiratory interaction. The anal- 
ysis of physiological signals of this kind has already been 
found useful in relation to several different diseases and phys- 
iological states (see e.g. [42] and references therein). Transi- 
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Figure 9: Analysis of signals from an analogue simulation of the 
system (fT4l) . (a) Phase portrait from the oscilloscope; (b) frequency 
cDi(t) from the external signal generator; (c) detected frequency 
u)2(t) of the second driven oscillator; (d) Fast Fourier Transform 
(FFT) of the detected frequency u)2(t). 
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Figure 10: Synchronization and time- varying parameters in the car- 
diorespiratory interaction, (a) Standard 2:N synchrogram. (b) Syn- 
chronization index for ratios 2:8 and 2:9 as indicated, (c) Time- 
evolution of the cardiac fh(t) and respiratory f r (t) frequency. 
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Figure 1 1 : (Color online) Coupling functions in the cardiorespiratory 
interaction calculated at different times, (a)-(c) Coupling function 
fa) from the first oscillator, and (d)-(f) ^2(^1, fa) from the 
second oscillator. The window time intervals were calculated at: t — 
725 s for (a) and (d); t = 1200 s for (b) and (e); and at t = 1250 s 
for (c) and (f). The window length was t w — 50 s. 



tions in cardiorespiratory synchronization have been studied 
in relation to anaesthesia i43h and sleep cycles 10]. It is also 
known that modulations and time- varying sources are present, 
and that these can affect the synchronization between biolog- 
ical oscillators l42l |45[ 14611 . For comprehensive and reliable 
analysis a technique is needed that is able, not only to identify 
the time- varying information, but which will allow evaluation 
of the interacting measures (e.g. synchronization and direc- 
tionality), based solely on the information inferred from the 
signals. We will show that our technique meets these criteria. 

We analyse cardiorespiratory measurements from human 
subject under anaesthesia. Their breathing rate was held con- 
stant, being determined by a respirator. For such systems the 
analytic model is unknown, in contrast to analogue and nu- 
merical examples, but the oscillatory nature of the signal is 
immediately evident. The instantaneous cardiac phase was 
estimated by synchrosqueezed wavelet decomposition i47ll of 
the ECG signal. Similarly, the respiratory phase was extracted 
from the respiration signal. The final phase time- series were 
reached after protophase-phase transformation fl2h . 

Application of the inferential technique reconstructs the 
phase parameters that govern the interacting dynamics. Fig. 
fTUtc) shows the time-evolution of the cardiac and respira- 
tion frequencies. It is evident that the constant pacing of the 
breathing is well-inferred, and that the instantaneous cardiac 
frequency, i.e. "heart rate variability", increases with time. 
The inferred parameters, and their correlations, is used to de- 
tect the occurrence of cardiorespiratory synchronization and 
the corresponding synchronization ratio. The synchronization 
evaluation 7 sync = s(c) G {0,1}, shown in Fig. [TOtb) re- 
veals the occurrence of transitions between the synchronized 
and non- synchronized states, and transitions between differ- 
ent synchronization ratios: from 2:8 (i.e. 1:4) at the begin- 
ning to 2:9 in the later intervals. Because the evaluation of 
the synchronization state is based on all of the given details 
about the phase dynamics, the proposed method not only de- 
tects the occurrence of transitions, but also describes their in- 
herent nature. The results for 7 svnc were consistent with the 



corresponding synchrogram shown in Fig.[T3a), but provided 
a clearer and less ambiguous indication of synchronization. 

The functional relationships that describe the cardiorespira- 
tory interactions are shown in Fig.[TTJ Evaluated for three dif- 
ferent time windows, the upper figures (a)-(c) show the cou- 
pling function q\ (fa , fa) from the cardiac oscillating activity, 
and the lower figures (d)-(f) show q2(fa,fa) from the respi- 
ration oscillator. The form of the functions is complex, and it 
changes qualitatively over time - cf. Fig. \TJla) with (b) and 
(c), or (d) with (e) and (f). The influence from respiration to 
heart (qi(fa,fa)) has a larger norm (i.e. coupling strength) 
than in the opposite direction, indicating that the predominant 
direction of coupling is from respiration to heart. One can also 
observe that q\ (fa ,fa) in (b) and (c) is of a fairly regular si- 
nusoidal form with a strong influence from respiration. This 
arises from the contribution of those base functions describing 
the direct influence of respiration (for a detailed discussion see 
Ugh). Furthermore, Fig.[TT]shows that the functional relation- 
ships for the interactions of an open (biological) system can 
in themselves be time- varying processes. 



VII. GENERALIZATION TO NETWORKS OF 
OSCILLATORS 

Our parameter inference procedure can be applied with 
only minimal modification to any number N of interacting 
oscillators within a general coupled-network structure. 

The notation of Eq. (Q]) is readily generalized for the N 
oscillators, and the inference procedure, Eq. ©, is then ap- 
plied to the corresponding iV-dimensional phase observable. 
The number of base function should of course be refined to 
include terms relevant for the iV-oscillator inference. Sub- 
procedures like the time-varying propagation, and the noise 
inference, apply exactly as before. We note that the computa- 
tional power required increases exponentially with N, which 
makes the method unsuitable for the inference of large-scale 
networks. However, for relatively small networks, a standard 
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Figure 12: (Color online) Coupling functions for the interacting network i\5l . Each row represents the influence on a specific oscillator: (a)-(c) 
on the first one, (d)-(f) on the second, and (g)-(i) on the third. The notation is such that e.g. ^23 represents the influence of the third oscillator 
on the second, while ^213 represents the join influences of the first and third oscillators on the second one. The numbers on the right of each 
coupling function represent their normalized coupling strength c and the significance p- value. The significant couplings are denoted with black 
squares. The parameter values were: uj\ = 2ir 1.1, 002 = 2ty 0.27, 003 = 2tt 3, £12 = 0.3, £13 = 0.2, £21 = 0.2, £312 = 0.5 and noise strength 
Ei — 0.1. The phases were estimated as fa = arctan(^/x;). 



high-performance computer will suffice for useful inference. 

We first demonstrate the inference on three interacting 
Poincare oscillators subject to noise 

3 jk 

V% -nyi +UiXi ^^eijyj + ^2e ijk yjyk + &(*), 

3 jk 

r i = (y/x 2 i +y?-l) i,j,k = 1,2,3, 

(15) 

where many of the coefficients and are initially set 
to zero; but some are non-zero, such as when the first oscil- 
lator is pairwise coupled to the second and third oscillators. 
The second oscillator is coupled also to the first (forming a 
bidirectional interaction). The third oscillator is influenced 
by the join contribution from the first and second oscillators. 
The latter coupling means physically that part of the network 
(cluster) exhibits a common functional influence on the other 
oscillators. The inference of this cross-coupling is the direct 
benefit of network (rather than pairwise) coupling detection. 

The inference of the three-dimensional phase variables 
from a numerical simulation of the network (IT31) is presented 
in Fig. [T2j The plots present the specific forms of coupling 
function that govern the interactions within the network. The 
coupling strengths are evaluated as partial norms from the rel- 
evant base functions. Note that the cross-couplings (c), (f) and 
(i) are shown for visual presentation as functions dependent on 



two phases, whereas the coupling strengths include also the 
base function dependent on all the three phases. In order to 
determine whether the inferred couplings are real or spurious, 
we conducted surrogate testing. The detected couplings were 
tested for significance in respect of 100 cou plin gs evaluated 
from surrogate phases. Cyclic surrogates |49|,[50|] were gener- 
ated from each of the phases, randomizing the temporal cross 
correlations, while preserving the frequencies and statistical 
characteristics unchanged. 

Recently, Kralemann et al. 115 ill discussed the notion of ef- 
fective and structural connectivity in networks. Effective cou- 
plings are those that are detected, while not present in the orig- 
inal structure e.g. indirectly-induced coupling. The question 
posed was: are the effective couplings real, or are they arti- 
facts? Our analysis showed that when one applies appropri- 
ate surrogate testing, the technique is able to distinguish the 
structural couplings as being significant. The resultant cou- 
pling strengths and significance p- values in Fig. [121 suggest 
that the connectivity (black-boxed couplings) of the network 
(IT31) was inferred correctly. Note also that some relations as 
in (g) have relatively large strength, even though they are less 
significant then some lower couplings as in (e). If the pos- 
sibility of effective couplings cannot be excluded, then our 
technique (with use of surrogate testing) provides a consistent 
way of inferring the true structure of the network. 

More importantly, the use of our method allows one to 
follow the time-variability of the structural and functional 
connectivity within the network. This is especially impor- 
tant when inferring the interactions of biological oscillators, 
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Figure 13: (Color online) Inference of time- varying coupling struc- 
ture for the network dT6t . The color/grayscale code for the couplings 
is presented in the box at the top. The four couplings £13, £14, £324 
and £42 were held constant at different values within three time seg- 
ments each of length 500 s. However, £21 was varied continuously 
through the whole time interval. For each segment the structure of 
the network is presented schematically on the diagrams in the dashed 
grey boxes. The parameters are given in the text. 

for which it is known that the dynamics is time-varying 
To illustrate the latter we infer the following 
network of four phase oscillators subject to white Gaussian 
noise 

0i =u 1 +asin(0i) + £13 (t) sin(0 3 ) + £i 4 (t) sin(0 4 ) 
(j) 2 = u 2 + asin(0 2 ) +e 2 i(t) sin(0 2 - <t>\) + 6(*) 

03 = ^3 + asin(0 3 ) + £324(2) sin(0 2 - 4 ) + 

04 = ^4 + asin(0 4 ) + £42^) sin(0 2 ) + U(t) • 

(16) 

Note that, because the coupling strengths are functions of 
time, we were effectively changing the structural connectivity 
of the network by varying their values. The parameter val- 
ues for the simulations were: uj\ = 27rl.ll, uj 2 = 27r2.13, 
ujs = 27r2.97, uj\ = 27r0.8, a = 0.2, and noise strengths 
Ei = 0.1. The couplings were varied discreetly in three time- 
segments, as follows, (i) For 0-500s: £13 = 0.4, £14 = 0.0, 
£324 = 0.4 and £42 = 0.4. (ii) For 500-1000s: £13 = 0, 
£14 = 0.35, £324 = and e 42 = 0.4. (iii) For 1000-1500: 
£13 = 0.45, £14 = 0.35, £324 = and £42 = 0. The coupling 
£21 was continuously varied between 0.5 — >• 0.3. The results 
presented in Fig. [T3l demonstrate that the method follows the 
time- variability of the couplings effectively and precisely. The 
dynamical variations are taking the network structure through 
various different connectivity states, and the different topolo- 
gies are detected reliably throughout their time-evolution. 

VIII. CONCLUSIONS 

Starting from the perspective of dynamical systems infer- 
ence, we have built an algorithm able to detect synchroniza- 



tion, to describe the functional form of the mutual interactions 
between oscillators, and to perform such tasks successfully in 
the presence of a time-evolving dynamics. 

The algorithm differs substantially from earlier approaches 
with respect both to synchronization detection capabilities, 
and to the estimation of coupling and directionality. Most 
other techniques are based on information flow (e.g. trans- 
fer entropy, or Granger causality) providing them with great 
generality. While limiting ourself to the hypothesis of con- 
tinuous time differential equations driving the dynamics (cor- 
respondingly restricting the domain of applicability), we can 
optimally exploit the benefits of this assumption. Unlike all 
other approaches, our technique does not require the observ- 
able to fill the domain of the probability density function at 
equilibrium. Thus, in oscillators with a limit cycle (Van der 
Pol, Fitz-Hugh Nagumo, Poicare, etc..) even one single ex- 
treme path is sufficient to characterize the parameters of the 
dynamics. Hence, we can determine uniquely the limit- time 
equilibrium distribution, i.e. the Fokker-Plank equation asso- 
ciated with the SDE. Thus an immediate advantage is that we 
can extract the same information from a fraction of the volume 
of data that is typically required by earlier methods. Because 
a very wide range of natural and artificial systems are describ- 
able in terms of continuous time differential equations (e.g. 
oscillatory processes in nature, mechanical systems, analogue 
voltage systems), the loss of generality in our approach is ac- 
tually minimal, compared to the advantage gained in terms of 
informational efficiency. 

We have applied the algorithm successfully to a representa- 
tive classes of oscillators, testing it on synthetically-generated 
data created from various models, and on data from an ana- 
logue circuit device with known dynamics. In each case, we 
were able to demonstrate the precision of parameter detection, 
the temporal precision of synchronization detection, and the 
accuracy of directionality identification. 

We have also demonstrated the efficacy of the technique in 
relation to cardiorespiratory time series data. Synchronization 
phenomena were already well-known in such systems, but the 
details of functional coupling were not. From the inferred pa- 
rameters we were able to reconstruct the extent of the cardiac 
and respiratory variability, estimate the direction of coupling, 
and detect the presence of and type of intermittent cardiores- 
piratory synchronization. 

Because the whole enterprize is built on an inference algo- 
rithm for an iV-dimensional dynamical system, the technique 
was readily extensible to the study of a network of oscillators 
whose parameters and coupling functions may be changing in 
time. An example of such application we considered a net- 
work of Poincare oscillators, generated by numerical simula- 
tion. We were able to demonstrate effective coupling detec- 
tion, cross- validating the results by surrogate testing. 

Although the implementation itself might see future im- 
provements (e.g. in terms of speed of calculation, or auto- 
matic base function selection), it is worth emphasizing that 
the method allows one to designate which components of the 
system are expected to be time- variable. Such selection is op- 
tional, but it provides an effective means by which to incor- 
porate previous knowledge available for any particular sys- 
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tern, and enables the algorithm to adapt itself optimally to the 
externally-imposed constraints. 

Given the advantages that the dynamical approach offers 
in tackling synchronization detection and coupling identifica- 
tion, we believe that the framework presented above will be 
found valuable for a wide range of future applications. 

Appendix A: Fixed point algorithmic check 

The procedure of synchronization detection between two 
oscillators generating phase time-series reduces to the investi- 
gation of synchronization of the synthetic phase model using 
the parameters returned by the Bayesian algorithm. To cal- 
culate s(c) for any of the sampled parameter sets, one can 
proceed as follows: 

(i) From an arbitrary fixed £, and for an arbitrary ^o, 
integrate numerically (using the standard fourth-order 
Runge-Kutta algorithm) the dynamical system pre- 
scribed by the phase base function (Eq. © without the 
noise) for one cycle of the toroidal coordinate, obtain- 



ing the mapped point M(ipo). 

(ii) Repeat the same integration for multiple ^ coordinates 
next to the initial one, obtaining the map M(ipi) 

(iii) Based on finite difference evaluation of dM/difj, use 
a modified version of Newton's root-finding method to 
analyse the function M(ip)—ip. The method is modified 
by calculating M at the next point ^n+i such that 

V>n+1 = + 0.8 X |(M(^n) - ^n)/(M'(^n) - 1))|- 

The coefficient 0.8 is an arbitrary constant that we 
found to be particularly efficient for solution of the 
problem. Note that in this version, Newton's method 
can only test the function by moving forward; in actual 
fact (a) the existence of the root is not guaranteed; and 
(b) we are not interested in the root itself but only in its 
existence. 

iv) If there is a root, s(c) = 1 is returned. If a root is not 
found, s(c) = is returned. 
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